{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: \"Anamak\" tokamak\n", "---\n", "\n", "Here, we generate a few different equilibria using the \"Anamak\" example tokamak from the [FGE code paper](https://arxiv.org/abs/2512.06847).\n", "\n", "The equilbirium\\profile parameters are **completely made up** based on some equilbiria in the above paper - please experiment on your own and change them to more realistic values as you please!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the machine object\n", "\n", "The Anamak machine is a useful toy tokamak because the geometry can be varied in a simple and controlled way. In particular, the number and placement of active coils, passive structures, and the vessel can all be adjusted using a small set of geometric parameters.\n", "\n", "Following the parameterisation used in the reference above, the $N_{\\mathrm{coils}}$ **active coils** are placed on an **ellipse** centred at $(R_0, 0)$ in the poloidal $(R, Z)$ plane. Their coordinates are defined by\n", "\n", "$$\n", "R_{\\mathrm{coil}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{coil}} \\cos(\\theta_i),\n", "$$\n", "\n", "$$\n", "Z_{\\mathrm{coil}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{coil}} \\sin(\\theta_i),\n", "$$\n", "\n", "where\n", "- $R_0$ is the major radius at which the ellipse is centred.\n", "- $\\hat{R}_{\\mathrm{coil}}$ is the minor radius of the coils. \n", "- $\\bold{\\theta}$ is a vector (length $N_{\\mathrm{coils}}$) of poloidal angles $\\theta \\in [0, 2\\pi)$. \n", "- $\\kappa$ is the \"elongation\" of the ellipse.\n", "\n", "Similarly, the $N_{\\mathrm{passives}}$ **passive structures** are placed on an ellipse at the same centre:\n", "\n", "$$\n", "R_{\\mathrm{passives}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{passive}} \\cos(\\theta_i),\n", "$$\n", "\n", "$$\n", "Z_{\\mathrm{passives}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{passive}} \\sin(\\theta_i),\n", "$$\n", "\n", "but with a potentially different minor radius $ \\hat{R}_{\\mathrm{passive}}$. \n", "\n", "Finally the $N_{\\mathrm{limiter}}$ limiter/wall points are set up similarly:\n", "\n", "$$\n", "R_{\\mathrm{limiter}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{limiter}} \\cos(\\theta_i),\n", "$$\n", "\n", "$$\n", "Z_{\\mathrm{limiter}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{limiter}} \\sin(\\theta_i).\n", "$$\n", "\n", "We should note that given we are only doing static solves below, the passive structures are not really used here." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell, we can choose the parameters for building the Anamak." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# major radius of the machine\n", "R0 = 1.0\n", "\n", "# \"elongation\" of the machine (same for actives, passives, and limiter here)\n", "kappa = 1.0 \n", "\n", "# for each of the coils, passives, and limiter, define:\n", "# -> N = the number of each you desire\n", "# -> R = the minor radius of its ellipse\n", "\n", "# actives parameters\n", "N_actives = 8\n", "R_coil = 0.7\n", "\n", "# passives parameters\n", "N_passives = 50\n", "R_passive = 0.5\n", "\n", "# limiter parameters\n", "N_limiter = 101\n", "R_limiter = 0.45" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following functions are used to build the Anamak." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# builds an ellipse of (R,Z) points\n", "def build_ellipse(\n", " R0, # major radius\n", " N, # number of points in [0,2\\pi)\n", " R_minor, # minor radius\n", " kappa, # \"elongation\"\n", " ):\n", " \n", " # poloidal locations\n", " theta = np.arange(0, 2*np.pi, 2*np.pi/N)\n", " \n", " # locations\n", " r_locs = R0 + R_minor*np.cos(theta)\n", " z_locs = R_minor*kappa*np.sin(theta)\n", "\n", " return r_locs, z_locs\n", "\n", "# builds active coils using (R,Z) locations of filaments\n", "def build_actives(\n", " r_locs, \n", " z_locs, \n", " ):\n", "\n", " # populate active coils dictionary with single filament coils (point-sources)\n", " active_coils_data = {}\n", " for i, (r, z) in enumerate(zip(r_locs, z_locs)):\n", " active_coils_data[i] = {\n", " \"R\": [r],\n", " \"Z\": [z],\n", " \"dR\": 0.05, # size only used for plotting in the simulations below\n", " \"dZ\": 0.05, # size only used for plotting in the simulations below\n", " \"resistivity\": 1.55e-8, # not used below\n", " \"polarity\": 1, # default\n", " \"multiplier\": 1, # single turn coils\n", " }\n", " \n", " return active_coils_data\n", "\n", "# builds passive structures using (R,Z) locations of filaments\n", "def build_passives(\n", " r_locs, \n", " z_locs, \n", " ):\n", "\n", " # populate passive coils list with single filament passives (point-sources)\n", " passive_coils_data = []\n", " for r, z in zip(r_locs, z_locs):\n", " passive_coils_data.append({\n", " \"R\": r,\n", " \"Z\": z,\n", " \"dR\": 0.02, # size only used for plotting in the simulations below\n", " \"dZ\": 0.02, # size only used for plotting in the simulations below\n", " \"resistivity\": 5.5e-7 # resistivity of material (steel)\n", " })\n", " \n", " return passive_coils_data\n", "\n", "# builds limiter geometry using (R,Z) locations\n", "def build_limiter(\n", " r_locs, \n", " z_locs, \n", " ):\n", "\n", " # populate limiter list\n", " limiter_data = []\n", " for r, z in zip(r_locs, z_locs):\n", " limiter_data.append({\"R\": r, \"Z\": z})\n", " \n", " return limiter_data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build the active coils data\n", "r_active_locs, z_active_locs = build_ellipse(\n", " R0=R0, \n", " N=N_actives, \n", " R_minor=1.4*R_passive,\n", " kappa=kappa,\n", ")\n", "\n", "active_coils_data = build_actives(\n", " r_locs=r_active_locs, \n", " z_locs=z_active_locs\n", " )\n", "\n", "# build the passive structures data\n", "r_passive_locs, z_passive_locs = build_ellipse(\n", " R0=R0, \n", " N=N_passives, \n", " R_minor=R_passive,\n", " kappa=kappa,\n", ")\n", "\n", "passive_coils_data = build_passives(\n", " r_locs=r_passive_locs, \n", " z_locs=z_passive_locs\n", " )\n", "\n", "# build the limiter data\n", "r_limiter_locs, z_limiter_locs = build_ellipse(\n", " R0=R0, \n", " N=N_limiter, \n", " R_minor=R_limiter,\n", " kappa=kappa,\n", ")\n", "\n", "limiter_data = build_limiter(\n", " r_locs=r_limiter_locs, \n", " z_locs=z_limiter_locs\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build the Anamak\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_data=active_coils_data,\n", " passive_coils_data=passive_coils_data,\n", " limiter_data=limiter_data,\n", " wall_data=limiter_data,\n", ")\n", "\n", "# plot the machine\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "tokamak.plot(axis=ax1,show=False) \n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) \n", "ax1.set_xlabel(\"Major radius [m]\")\n", "ax1.set_ylabel(\"Height [m]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate an equilibrium" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import equilibrium_update\n", "\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=1.9, # radial range\n", " Zmin=-0.8, Zmax=0.8, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate a profile object\n", "\n", "We use the Lao profile object here, providing constant pressure and toroidal field profiles for simplicity. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke.jtor_update import Lao85\n", "profiles = Lao85(\n", " eq=eq, # eq object\n", " Ip=2e5, # plasma current\n", " fvac=0.5*R0, # toroidal field * major radius\n", " alpha=[58213.6], # linear profile\n", " beta=[0.582136] # linear profile\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the static nonlinear solver" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Inverse solver constraints\n", "\n", "In this first example, we'll set up some simple constraints to generate an up/down symmetric double null plasma (with relatively high triangularity). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# null points constraints\n", "Rx = 0.8 # X-point radius\n", "Zx = 0.35 # X-point height\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# isoflux constraints (including null points here)\n", "Rout = 1.22 # outboard midplane radius\n", "isoflux_set = np.array([\n", " [\n", " [Rx, Rx, Rout], \n", " [Zx, -Zx, 0.0]\n", " ]]\n", " )\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", " weight_isoflux=1.0,\n", " weight_nulls=0.8, # slightly lower weight on null points to get correct shape\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The inverse solve\n", "\n", "Solve using the inverse solver with forced up/down symmetry.\n", "\n", "Setting `force_up_down_symmetric=True` projects the initial plasma flux onto its even-in-Z component and keeps the plasma response even in every forward solve, including those used to construct the full inverse Jacobian. The grid, machine, and requested constraints should also be up/down symmetric. This flag constrains the plasma response, not the coil-current optimisation: if exactly equal currents are required in separate upper and lower coils, those currents must additionally be linked or constrained by the machine/control model." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-5,\n", " target_relative_psit_update=1e-3,\n", " verbose=True,\n", " l2_reg=np.array([1e-12]*N_actives),\n", " force_up_down_symmetric=True,\n", " )\n", "\n", "plasma_flux_parity_error = np.linalg.norm(eq.plasma_psi - eq.plasma_psi[:, ::-1]) / np.linalg.norm(eq.plasma_psi)\n", "print(f\"Relative plasma-flux parity error = {plasma_flux_parity_error:.3e}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# refine with a forward solve (now the currents are known)\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9,\n", " verbose=False, # print output\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using only a few constraints, we can see that the inverse solver has found the desired plasma. It is almost perfectly up/down symmetric and has high triangularity. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Triangularity = {eq.triangularity()}\")\n", "print(f\"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll use the same machine description to make a negtive triangularity plasma (again up/down symmetric). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# null points constraints\n", "Rx = 1.1 # X-point radius\n", "Zx = 0.35 # X-point height\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# isoflux constraints (including null points here)\n", "Rin = 0.8 # inboard midplane radius\n", "isoflux_set = np.array([\n", " [\n", " [Rx, Rx, Rin], \n", " [Zx, -Zx, 0.0]\n", " ]]\n", " )\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", " weight_isoflux=1.0,\n", " weight_nulls=0.8, # slightly lower weight on null points to get correct shape\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solve\n", "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-5,\n", " target_relative_psit_update=1e-3,\n", " verbose=True,\n", " l2_reg=np.array([1e-14]*N_actives),\n", " force_up_down_symmetric=True,\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# refine with a forward solve (now the currents are known)\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9,\n", " verbose=False, # print output\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Triangularity = {eq.triangularity()}\")\n", "print(f\"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets look at changing the shape of the Anamak.\n", "\n", "We'll now increase the number of active coils and increase the \"elongation\" of the device. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# \"elongation\" of the machine (same for actives, passives, and limiter here)\n", "kappa = 1.8\n", "\n", "# for each of the coils, passives, and limiter, define:\n", "# -> N = the number of each you desire\n", "# -> R = the minor radius of its ellipse\n", "\n", "# actives parameters\n", "N_actives = 24\n", "R_coil = 0.7" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build the active coils data\n", "r_active_locs, z_active_locs = build_ellipse(\n", " R0=R0, \n", " N=N_actives, \n", " R_minor=1.4*R_passive,\n", " kappa=kappa,\n", ")\n", "\n", "active_coils_data = build_actives(\n", " r_locs=r_active_locs, \n", " z_locs=z_active_locs\n", " )\n", "\n", "# build the passive structures data\n", "r_passive_locs, z_passive_locs = build_ellipse(\n", " R0=R0, \n", " N=N_passives, \n", " R_minor=R_passive,\n", " kappa=kappa,\n", ")\n", "\n", "passive_coils_data = build_passives(\n", " r_locs=r_passive_locs, \n", " z_locs=z_passive_locs\n", " )\n", "\n", "# build the limiter data\n", "r_limiter_locs, z_limiter_locs = build_ellipse(\n", " R0=R0, \n", " N=N_limiter, \n", " R_minor=R_limiter,\n", " kappa=kappa,\n", ")\n", "\n", "limiter_data = build_limiter(\n", " r_locs=r_limiter_locs, \n", " z_locs=z_limiter_locs\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build the Anamak\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_data=active_coils_data,\n", " passive_coils_data=passive_coils_data,\n", " limiter_data=limiter_data,\n", " wall_data=limiter_data,\n", ")\n", "\n", "# plot the machine\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "tokamak.plot(axis=ax1,show=False) \n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) \n", "ax1.set_xlabel(\"Major radius [m]\")\n", "ax1.set_ylabel(\"Height [m]\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Re-initialise the equilibrium and profiles objects: here we've increased the vertical size of the computational grid and increased the plasma current. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import equilibrium_update\n", "\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=0.1, Rmax=1.9, # radial range\n", " Zmin=-1.5, Zmax=1.5, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", ")\n", "\n", "from freegsnke.jtor_update import Lao85\n", "profiles = Lao85(\n", " eq=eq, # eq object\n", " Ip=5e5, # plasma current\n", " fvac=0.5*R0, # toroidal field * major radius\n", " alpha=[58213.6], # linear profile\n", " beta=[0.582136] # linear profile\n", ")\n", "\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try to create a simple up/down symmetric double null plasma in this new machine. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# null points constraints\n", "Rx = 0.8 # X-point radius\n", "Zx = 0.65 # X-point height\n", "null_points = [[Rx, Rx], [Zx, -Zx]]\n", "\n", "# isoflux constraints (including null points here)\n", "Rout = 1.3 # outboard midplane radius\n", "isoflux_set = np.array([\n", " [\n", " [Rx, Rx, Rout], \n", " [Zx, -Zx, 0.0]\n", " ]]\n", " )\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", " weight_isoflux=1.0,\n", " weight_nulls=0.7, # slightly lower weight on null points to get correct shape\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solve\n", "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-5,\n", " target_relative_psit_update=1e-3,\n", " verbose=True,\n", " l2_reg=np.array([1e-12]*N_actives),\n", " force_up_down_symmetric=True,\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# refine with a forward solve (now the currents are known)\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9,\n", " verbose=False, # print output\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us try to use the snowflake (i.e. second order null) constraints to make something a bit more exotic. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# null points constraints\n", "Rx = 0.8 # X-point radius\n", "Zx = 0.65 # X-point height\n", "null_points_2nd_order = [[Rx, Rx], [Zx, -Zx]] \n", "\n", "# isoflux constraints (including null points here, not really needed but should be fine)\n", "isoflux_set = np.array([\n", " [\n", " [Rx, Rx], \n", " [Zx, -Zx]\n", " ]]\n", " )\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(\n", " null_points_2nd_order=null_points_2nd_order,\n", " # isoflux_set=isoflux_set,\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solve\n", "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-5,\n", " target_relative_psit_update=1e-3,\n", " verbose=True,\n", " Picard_handover=1e-5,\n", " l2_reg=np.array([1e-12]*N_actives),\n", " force_up_down_symmetric=True,\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now perhaps a single null plasma - note we change some solver settings here too. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# null points constraints\n", "Rx = 0.8 # X-point radius\n", "Zx = 0.6 # X-point height\n", "null_points = [[Rx], [-Zx]]\n", "\n", "# isoflux constraints (including null points here)\n", "Rout = 1.32 # outboard midplane radius\n", "isoflux_set = np.array([\n", " [\n", " [Rx, Rout, 0.85, 1.17, 0.65], \n", " [-Zx, 0.0, 0.57, 0.45, 0.0]\n", " ]]\n", " )\n", " \n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", " weight_isoflux=1.0,\n", " weight_nulls=0.7, # slightly lower weight on null points to get correct shape\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solve\n", "eq.set_plasma_psi(eq.create_psi_plasma_default())\n", "eq.tokamak.set_all_coil_currents(current_vec=np.zeros(N_actives))\n", "\n", "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-5,\n", " target_relative_psit_update=1e-3,\n", " verbose=True,\n", " l2_reg=np.array([1e-14]*N_actives),\n", " Picard_handover=1e-5,\n", " full_jacobian_handover=[2e-2,1e-2],\n", " force_up_down_symmetric=False,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints" ] } ], "metadata": { "kernelspec": { "display_name": "freegsnke", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.20" } }, "nbformat": 4, "nbformat_minor": 4 }